A good way to introduce OOP is through a simple example that uses Microsoft Excel
Here, we see two main components: a) some data (the values in the spreadsheet), and b) a request to compute an average (applied to the data in column B, rows 2 through 6)
One can think of of these components more generally as, respectively, a) an object, and b) a function
An object is a facility that stores data or information, while a function is a prescribed set of actions that will be applied to one or more aspects of an object
Note that a function can also be considered an object
R and other OOPs abstract from environments like those seen in Excel (or Stata, etc.), allowing users to flexibly construct, manipulate, and deploy many, many objects and functions at the same time
Learning how to use R is all about learning how to work with different types of R objects and functions, and then eventually using these lessons to further your own research and statistical needs
The R and RStudio Environments
Essential Elements:
Console (R/RStudio)
Scripts (R/RStudio)
Environment/History (RStudio)
Files/Plots/Packages/Help (RStudio)
Tab Completion
Highlight and In-line Execution
Getting Help
Internal R methods include help() or ?, as well as ??
R can act just like a scientific computer, via very predictable commands
These commands can be roughly divided between those that are arithmetic (yielding numerical results) and those that are boolean (yielding logical results)
Though R can carry out the same operations as a calculator, what sets it apart is its ability to store calucated results (or any data or information) as objects
Three main assignment operators are used to populate objects:
<-
->
=
Conventional use is of the format objectname <- objectdata
mysum <-4 +5
simple.division <-4/5
Several semi-complex naming conventions apply to object names, but you generally won't have a problem if your object names a) start with letters and b) only use . and _ from the punctuation marks
Numbers are fine, so long as they do not begin the object name
15myobject<-4+5
myobject<-4+5
The vast majority of punctuation marks are reserved for other operations (~, !, @, #, $, %, &, *, ;, ', ?, |)
$myobject<-4+5
myobject<-4+5
Note also that object names and data in R are case-sensitive
myobject <-4 +5print(MYOBJECT)
print(myobject)
Assignment Operators/Concatenation
Often we want to assign multiple value to an object
To do this, we need to use the concatenation operator, or c()
mynumbers <-c(4, 5, 6)
print(mynumbers)
## [1] 4 5 6
We can also take advantage of some shortcuts for assigning a large number of values to an object, or when our values of interest following a particular pattern
In OOP languages like R, objects that take a certain form or have particular attributes can be assigned to a class
The class of an object acts like a shortcut for many functions, allowing them to carry out certain calculations that are only appropriate for objects of a certain kind
The most prominent object classes in R include:
Numeric
Character
Factor
Logical
Missing/Non-numeric/Null
x <-4 +5class(x)
is.numeric(x)
y <- "New York"class(y)
is.character(y)
is.numeric(y)
z <-factor(x =c("New York", "New Jersey"), ordered = T)
print(z)
class(z)
is.factor(z)
is.character(z)
a <-TRUEis.logical(a)
is.character(a)
is.factor(a)
b <-NAclass(b)
is.na(b)
is.numeric(b)
c <-0/0class(c)
print(c)
is.nan(c)
is.na(c)
d <-5/0class(d)
is.infinite(d)
e <-NULLclass(e)
is.null(e)
is.numeric(e)
Object Classes/Data Structures
There are six basic data structures in R, or classes of objects that are commonly used to store data, namely:
Atomic Vectors -- a single piece of data
x <-3print(x)
y <- "New York"print(y)
Vectors -- a uni-dimensional ordered collection of class-homogenous data
x <-c(3, 4, 9, 10)
print(x)
y <-c("New York", "New Jersey", "New Mexico")
print(y)
z <-c(1, 2, "a", "b")
print(z)
class(z)
as.numeric(z)
Matrices -- a two-dimensional ordered collection of class-homogenous data
x <-matrix(data =1:100, nrow =10)
print(x)
y <-matrix(data =c(1, 2, "a", "b"), nrow =2)
print(y)
Arrays -- a generalized n-dimensional ordered collection of class-homogenous data
x <-array(data =1:300, dim =c(10, 10, 3))
print(x)
Data Frames -- a two-dimensional ordered collection of class-heterogenous data
x <-data.frame(c(1, 2), c("a", "b"))
print(x)
x <-data.frame(col1 =c(1, 2), col2 =c("a", "b"))
print(x)
class(x)
class(x[, 1])
class(x[, 2])
Lists -- a generalized n-dimensional ordered collection of class-heterogenous data
Often we're not interested in looking at the entirety of an object, but rather at specific attributes or subsets of that object
All individual pieces of data in an R object are assigned an index, whose dimensionality is object specific (e.g. a cell in a matrix will have a two-component index, or a data point in a six-dimensional array will have a six-component index)
This index can be used to isolate aspects of R objects via brackets
Instead of inspecting aspects of an object based on their index, we often like to inspect only those aspects of a data object that satisfy certain conditions
R allows us to do this by constructing a boolean object (containing only TRUE or FALSE for each element of an object) and only returning those object elements whose boolean positional analog is TRUE
missing.letters1 <-function(text) {
data.split <-strsplit(x = text, split ="") # Split the text into letters
data.split <-unlist(data.split) # Convert the above result from a list to a vector
data.split <-tolower(data.split) # Convert all letter to lowercase
data.split <-unique(data.split) # Keep only the unique letters
data.split <-data.split[data.split != " "] # Get rid of spaces
missing <-letters[!letters %in%data.split] # Find the letters that are NOT in the text return(missing) # Return the missing letters
}
missing.letters1(text ="All we have to fear is fear itself")
## [1] "b" "c" "d" "g" "j" "k" "m" "n" "p" "q" "u" "x" "y" "z"
## Error: multi-argument returns are not permitted
In theory, R functions don't permit multi-argument returns. What we need to do is "trick" R into returning multiple arguments by handing it a sort of Matryoshka doll. The easiest way to do this is to put all our desired return arguments into a list, and then have our function return the list.
# Get the second element of the list
example.list[[2]]
## [1] 15
Programming/Functions/Conditional Execution
An essential element of programming in any language is the use of conditional execution statements, more commonly known as if/else statements. The syntax for these in R is pretty straight forward, taking the basic form of:
if (condition) true.action else false.action
Some examples with vectors of length == 1:
x <-10
if (x <10) "fine" else "also fine"
## [1] "also fine"
if (x >=10) x *100 else 50/x
## [1] 1000
if (x <10) {
y <-2 *x
c("really great", y)
} else {
y <--2 *x
c("really not great", y)
}
## [1] "really not great" "-20"
Things fall apart, however, when we deal with vectors of length > 1:
x <-8:12
if (x <10) "fine" else "also fine"
## Warning: the condition has length > 1 and only the first element will be
## used
## [1] "fine"
if (x >=10) x *100 else 50/x
## Warning: the condition has length > 1 and only the first element will be
## used
## [1] 6.250 5.556 5.000 4.545 4.167
if (x <10) {
y <-2 *x
c("really great", y)
} else {
y <--2 *x
c("really not great", y)
}
## Warning: the condition has length > 1 and only the first element will be
## used
For vectors of length > 1, we need to be sure to vectorize our conditional execution statements (meaning, they must be instructed to evaluate objects element-by-element). R has three basic facilities for doing this: ifelse(), for loops, and the family of apply functions
x <-1:1e+05
time.ifelse <-system.time(ifelse(test = x <10, yes = {
y <-2 *x
3 *y
}, no = {
y <--2 *x
3 *y
}))
out <-NULL
time.forloop <-system.time(for (i in 1:length(x)) {
cat(i, "\n")
if (x[i] <10) {
y <-2 *x[i]
out[i] <-3 *y
} else {
y <--2 *x[i]
out[i] <-3 *y
}
})
time.forloop/time.ifelse
Useful Stuff: Applied Research Edition
For most applied researchers, "useful stuff" that can be done in R boils down to a few core items:
Importing different types of data from different sources
Summarizing the structure and content of data
Carrying out operations and calculations across groups
Reshaping data to and from various formats
Attempting to conduct causal inference
Importing data/spreadsheets
For spreadsheet data that is not explicitly saved as a Microsoft Excel file:
# Import data with header row, values separated by ',', decimals as '.'
data <-read.csv(file =" ", stringsAsFactors = ) # Basic Call
# Example
cpds.data <-read.csv(file ="datasets/cpds_comma.csv", stringsAsFactors = F)
dim(cpds.data) # Check dimensionality of dataset
## [1] 1173 237
cpds.data[1:5, 1:5] # View first 5 rows/columns
## year country countryn gov_right gov_cent
## 1 1960 Australia 1 100 0
## 2 1961 Australia 1 100 0
## 3 1962 Australia 1 100 0
## 4 1963 Australia 1 100 0
## 5 1964 Australia 1 100 0
# Import data with header row, values separated by ';', decimals as ','
data <-read.csv2(file =" ", stringsAsFactors = )
# Import data with header row, values separated by tab, decimals as '.'
data <-read.delim(file =" ", stringsAsFactors = )
# Import data with header row, values separated by tab, decimals as ';'
data <-read.delim2(file =" ", stringsAsFactors = )
# For importing tabular data with maximum customizeability
data <-read.table(file = , header = , sep = , quote = , dec = , fill = , stringsAsFactors = )
For spreadsheet data that is explicitly saved as a Microsoft Excel file (.xls or .xlsx):
# Install the 'gdata' package (only necessary one time)# install.packages('gdata') # Not Run# Load the 'gdata' package (necessary every new R session)library(gdata)
# For importing both .xls and .xlsx files
data <-read.xls(xls = , sheet = ) # Basic Call
# Example with .xls (single sheet)
cpds.data <-read.xls(xls ="datasets/cpds_excel_old.xls")
cpds.data[1:5, 1:5]
## year country countryn gov_right gov_cent
## 1 1960 Australia 1 100 0
## 2 1961 Australia 1 100 0
## 3 1962 Australia 1 100 0
## 4 1963 Australia 1 100 0
## 5 1964 Australia 1 100 0
# Example with .xlsx (single sheet)
cpds.data <-read.xls(xls ="datasets/cpds_excel_new.xlsx")
cpds.data[1:5, 1:5]
## year country countryn gov_right gov_cent
## 1 1960 Australia 1 100 0
## 2 1961 Australia 1 100 0
## 3 1962 Australia 1 100 0
## 4 1963 Australia 1 100 0
## 5 1964 Australia 1 100 0
# Install the 'foreign' package (only necessary one time)# install.packages('foreign') # Not Run# Load the 'foreign' package (necessary every new R session)library(foreign)
read.xxxx(file = ) # Basic call, where .xxxx is the imported file's extension/client# For importing Stata files
data <-read.dta(file = )
# Example
cpds.data <-read.dta(file ="datasets/cpds_stata.dta")
cpds.data[1:5, 1:5]
## year country countryn gov_right gov_cent
## 1 1960 Australia 1 100 0
## 2 1961 Australia 1 100 0
## 3 1962 Australia 1 100 0
## 4 1963 Australia 1 100 0
## 5 1964 Australia 1 100 0
# For importing SPSS files
data <-read.spss(file = )
# For importing SAS files
data <-read.ssd(file = )
data <-read.xport(file = )
# For importing Fortran files
data <-read.fortran(file = )
# For importing DBF files
data <-read.dbf(file = )
Importing data/urls (e.g.: http, https, ftp)
Most data importing facilities in R can be adapted to import non-local files via http/s/ftp
Importing data/other (e.g.: html tables, html, xml, json)
Say we wanted to scrape a list of votes in the US Congress contained in an html table (like this)
# Install the 'XML' package (only necessary one time)# install.packages('XML') # Not Run# Load the 'XML' package (necessary every new R session)library(XML)
data <-readHTMLTable(doc = , header = , which = , stringsAsFactors = ) # Basic Call
link <- "http://clerk.house.gov/evs/2014/ROLL_000.asp"
votes.2014 <-readHTMLTable(doc = link, header = T, which =1, stringsAsFactors = F)
dim(votes.2014) # Check data dimensionality
## Roll Date Issue
## 1 99 6-Mar H RES 501
## 2 98 5-Mar H R 2126
## 3 97 5-Mar H R 4118
## 4 96 5-Mar H R 4118
## 5 95 5-Mar H R 938
## 6 94 5-Mar H RES 497
## 7 93 5-Mar H RES 497
## 8 92 4-Mar H RES 488
## 9 91 4-Mar H R 3370
## 10 90 28-Feb H R 899
## 11 89 28-Feb H R 899
## 12 88 28-Feb H R 899
## 13 87 28-Feb H R 899
## 14 86 28-Feb H R 899
## 15 85 27-Feb H R 3193
## 16 84 27-Feb H R 3193
## 17 83 27-Feb H R 3193
## 18 82 27-Feb H R 3193
## 19 81 27-Feb H R 3193
## 20 80 27-Feb H RES 492
## 21 79 27-Feb H RES 492
## 22 78 27-Feb H R 2804
## 23 77 27-Feb H R 2804
## 24 76 27-Feb H R 2804
## 25 75 27-Feb H R 2804
## Question Result
## 1 On Ordering the Previous Question P
## 2 On Motion to Suspend the Rules and Pass, as Amended P
## 3 On Passage P
## 4 On Motion to Recommit with Instructions F
## 5 On Motion to Suspend the Rules and Pass, as Amended P
## 6 On Agreeing to the Resolution P
## 7 On Ordering the Previous Question P
## 8 On Motion to Suspend the Rules and Agree, as Amended P
## 9 On Motion to Suspend the Rules and Pass, as Amended P
## 10 On Passage P
## 11 On Motion to Recommit with Instructions F
## 12 On Agreeing to the Amendment F
## 13 On Agreeing to the Amendment F
## 14 On Agreeing to the Amendment F
## 15 On Passage P
## 16 On Motion to Recommit with Instructions F
## 17 On Agreeing to the Amendment F
## 18 On Agreeing to the Amendment A
## 19 On Agreeing to the Amendment A
## 20 On Agreeing to the Resolution P
## 21 On Ordering the Previous Question P
## 22 On Passage P
## 23 On Motion to Recommit with Instructions F
## 24 On Agreeing to the Amendment F
## 25 On Agreeing to the Amendment F
Or say we wanted to scrape the content of a web page (for example, here)
# Install the 'RCurl' package (only necessary one time)# install.packages('RCurl') # Not Run# Load the 'RCurl' and 'XML' packages (necessary every new R session)library(RCurl)
library(XML)
link <- "http://www.whitehouse.gov/the_press_office/President_Barack_Obamas_Inaugural_Address"
speech <-getURL(url = link)
# Optional html parsing (clean file, extract content)
speech.parsed <-htmlParse(file = speech)
text <-xpathSApply(doc = speech.parsed, path ="//div[@id='content']", fun = xmlValue)
# Optional content analysis# Install the 'qdap' package (only necessary one time)# install.packages('qdap') # Not Run# Load the 'qdap' package (necessary every new R session)library(qdap)
# Convert the text to a data.frame
text.df <-data.frame(text = text, stringsAsFactors = F)
# Split the text into sentences
sentences <-sentSplit(dataframe = text.df, text.var ="text")
# Compute various polarity (sentiment) statistics
polarity.score <-polarity(text.var = sentences$text)
polarity.score$group[, c("total.sentences", "total.words", "ave.polarity")]
Or say we wanted to scrape the content of an RSS feed (XML) (for example, here)
# Specify a link and issue a GET request
link <- "http://rss.nytimes.com/services/xml/rss/nyt/HomePage.xml"
page <-getURL(url = link)
# Parse the results of the get request as an XML file
rss.feed <-xmlParse(file = page) # Not printed
Or say we wanted to import a JSON file
# Install the 'jsonlite' package (only necessary one time)# install.packages('jsonlite') # Not Run# Load the 'jsonlite' package (necessary every new R session)library(jsonlite)
# Import JSON file, convert to list
cpds.data <-fromJSON(txt ="datasets/cpds_json.txt", simplifyDataFrame = F)
head(cpds.data)
## year country variable value
## 1 1960 Australia gov_right 100
## 2 1961 Australia gov_right 100
## 3 1962 Australia gov_right 100
## 4 1963 Australia gov_right 100
## 5 1964 Australia gov_right 100
## 6 1965 Australia gov_right 100
Summarizing data/example dataset
Replication data from Gelman et al., "Rich State, Poor State, Red State, Blue State: What’s the Matter with Connecticut?", Quarterly Journal of Political Science, 2007, 2: 345-367
For decades, the Democrats have been viewed as the party of the poor, with the Republicans representing the rich. Recent presidential elections, however, have shown a reverse pattern, with Democrats performing well in the richer blue states in the northeast and coasts, and Republicans dominating in the red states in the middle of the country and the south. Through multilevel modeling of individual- level survey data and county- and state-level demographic and electoral data, we reconcile these patterns.Furthermore, we find that income matters more in red America than in blue America. In poor states, rich people are much more likely than poor people to vote for the Republican presidential candidate, but in rich states (such as Connecticut), income has a very low correlation with vote preference.
# View last six rows and all columnstail(red.blue)
## state pres04 sex race age9 partyid income
## 76200 50 2 female white 65-74 republican under $15,000
## 76201 50 2 female white 50-59 independent under $15,000
## 76202 50 2 male white 18-24 independent under $15,000
## 76203 50 3 female black 30-39 something else under $15,000
## 76204 50 2 female white 50-59 democrat $30,000-$49,999
## 76205 50 2 male white 45-49 republican $15,000-$29,999
## relign8 age60 age65 geocode sizeplac
## 76200 catholic 60 or over 65 or over 1 city: 10,000 to 49,999
## 76201 protestant 45-59 50-64 1 rural
## 76202 other christian 18-29 18-24 1 city: 10,000 to 49,999
## 76203 other christian 30-44 30-39 1 city: 10,000 to 49,999
## 76204 mormon/lds 45-59 50-64 1 city: 10,000 to 49,999
## 76205 something else 45-59 40-49 2 rural
## brnagain attend year region y
## 76200 <NA> <NA> 2004 4 1
## 76201 <NA> <NA> 2004 4 1
## 76202 <NA> <NA> 2004 4 1
## 76203 <NA> <NA> 2004 4 NA
## 76204 <NA> <NA> 2004 4 1
## 76205 <NA> <NA> 2004 4 1
# Get detailed column-by-column informationstr(red.blue)
## 'data.frame': 76205 obs. of 17 variables:
## $ state : int 2 2 2 2 2 2 2 2 2 2 ...
## $ pres04 : int 1 2 1 1 1 1 1 2 2 2 ...
## $ sex : Factor w/ 2 levels "male","female": 2 1 2 2 2 2 1 2 2 2 ...
## $ race : Factor w/ 5 levels "white","black",..: 1 1 2 2 1 1 1 1 1 1 ...
## $ age9 : Factor w/ 9 levels "18-24","25-29",..: 2 1 3 3 4 3 4 1 2 1 ...
## $ partyid : Factor w/ 4 levels "democrat","republican",..: NA NA NA NA NA NA NA NA NA NA ...
## $ income : Factor w/ 8 levels "under $15,000",..: NA NA NA NA NA NA NA NA NA NA ...
## $ relign8 : Factor w/ 8 levels "protestant","catholic",..: NA NA NA NA NA NA NA NA NA NA ...
## $ age60 : Factor w/ 4 levels "18-29","30-44",..: 1 1 2 2 2 2 2 1 1 1 ...
## $ age65 : Factor w/ 6 levels "18-24","25-29",..: 2 1 3 3 4 3 4 1 2 1 ...
## $ geocode : int 3 3 3 3 3 3 3 3 3 3 ...
## $ sizeplac: Factor w/ 5 levels "city over 500,000",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ brnagain: Factor w/ 2 levels "yes","no": NA NA NA NA NA NA NA NA NA NA ...
## $ attend : Factor w/ 6 levels "more than once a week",..: NA NA NA NA NA NA NA NA NA NA ...
## $ year : num 2004 2004 2004 2004 2004 ...
## $ region : num 4 4 4 4 4 4 4 4 4 4 ...
## $ y : num 0 1 0 0 0 0 0 1 1 1 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr " 6 Jun 2007 14:53"
## - attr(*, "formats")= chr "%8.0g" "%8.0g" "%8.0g" "%8.0g" ...
## - attr(*, "types")= int 251 251 251 251 251 251 251 251 251 251 ...
## - attr(*, "val.labels")= chr "stanum" "presak04" "sex" "race" ...
## - attr(*, "var.labels")= chr "state id" "in today's election for president, did you just vote for:" "are you:" "are you:" ...
## - attr(*, "version")= int 8
## - attr(*, "label.table")=List of 14
## ..$ stanum : Named int 2
## .. ..- attr(*, "names")= chr "alaska"
## ..$ presak04: Named int 0 1 2 3 9
## .. ..- attr(*, "names")= chr "did not vote" "kerry" "bush" "nader" ...
## ..$ sex : Named int 1 2
## .. ..- attr(*, "names")= chr "male" "female"
## ..$ race : Named int 1 2 3 4 5
## .. ..- attr(*, "names")= chr "white" "black" "hispanic/latino" "asian" ...
## ..$ age9 : Named int 1 2 3 4 5 6 7 8 9
## .. ..- attr(*, "names")= chr "18-24" "25-29" "30-39" "40-44" ...
## ..$ partyid : Named int 1 2 3 4
## .. ..- attr(*, "names")= chr "democrat" "republican" "independent" "something else"
## ..$ income : Named int 1 2 3 4 5 6 7 8
## .. ..- attr(*, "names")= chr "under $15,000" "$15,000-$29,999" "$30,000-$49,999" "$50,000-$74,999" ...
## ..$ relign8 : Named int 1 2 3 4 5 6 7 8
## .. ..- attr(*, "names")= chr "protestant" "catholic" "mormon/lds" "other christian" ...
## ..$ age60 : Named int 1 2 3 4
## .. ..- attr(*, "names")= chr "18-29" "30-44" "45-59" "60 or over"
## ..$ age65 : Named int 1 2 3 4 5 6
## .. ..- attr(*, "names")= chr "18-24" "25-29" "30-39" "40-49" ...
## ..$ geocode : Named int 1 2 3
## .. ..- attr(*, "names")= chr "juneau/fairbanks/rural" "anchorage" "anchorage-fairbanks corridor"
## ..$ sizeplac: Named int 1 2 3 4 5
## .. ..- attr(*, "names")= chr "city over 500,000" "city: 50,000 to 500,000" "suburbs" "city: 10,000 to 49,999" ...
## ..$ brnagain: Named int 1 2
## .. ..- attr(*, "names")= chr "yes" "no"
## ..$ attend : Named int 1 2 3 4 5 9
## .. ..- attr(*, "names")= chr "more than once a week" "once a week" "a few times a month" "a few times a year" ...
Summarizing data/object content
A good place to start with our data is to calculate summary statistics
Caution: Although we have some inherently ratio (continuous) variables, the authors binned the data to create ordinal variables (e.g.: Age and Income). Bonus: the bins have differing widths.
Take a look at the variable "age9"
# View unique values observed in 'age9' (note the factor levels)unique(x = red.blue$age9)
## [1] 25-29 18-24 30-39 40-44 45-49 65-74
## [7] 75 or over 50-59 <NA> 60-64
## Levels: 18-24 25-29 30-39 40-44 45-49 50-59 60-64 65-74 75 or over
Binned Age
True Age
1
18-24
2
25-29
3
30-39
4
40-44
5
45-49
6
50-59
7
60-64
8
65-74
9
75+
Some notes on computing summary statistics:
Note that these functions are sensitive to missing values (NA); you should be sure to specify na.rm=T to avoid errors
# Sample 100 times from the standard normal distribution
sample.data <-rnorm(n =100, mean =0, sd =1)
# Add some missing values to the sample
sample.data[c(1, 4, 16, 64)] <-NA# Attempt to calculate the sample mean (presence of NAs)mean(x = sample.data)
## [1] NA
# Remove missing values from the sample
sample.data <-sample.data[!is.na(sample.data)]
# Attempt to calculate the sample mean (absence of NAs)mean(x = sample.data)
## [1] -0.1291
These functions are also sensitive to the presence of factor variables; remove the factor levels to avoid errors (usually use one of as.vector(), as.character(), or as.numeric())
# Check for factor levels in 'age9'is.factor(x = red.blue$age9)
## [1] TRUE
# Attempt to calculate the sample mean ('age9' is factor)mean(x = red.blue$age9, na.rm = T)
## Warning: argument is not numeric or logical: returning NA
## [1] NA
# Remove factor levels in 'age9'
red.blue$age9 <-as.numeric(x = red.blue$age9)
# Check that there are no more factor levels in 'age9'is.factor(x = red.blue$age9)
## [1] FALSE
# Attempt to calculate the sample mean ('age9' is not factor)mean(x = red.blue$age9, na.rm = T)
# Subset the gender vector to include only males
males <-gender[gender == "male"]
head(males)
## [1] "male" "male" "male" "male" "male" "male"
# Compute the percentage of males using vector lengthslength(males)
## [1] 34236
length(gender)
## [1] 75868
length(males)/length(gender) *100
## [1] 45.13
This can get really complicated really quickly
Suppose we wanted to know how voting behavior in the 2004 Presidential Election varies by race
That means we have to calculate frequencies as above for each race in {White, Black, Hispanic/Latino, Asian, Other} and each vote choice in {Bush, Kerry, Nader, Other, No Vote}.
How to tackle these tabulations?
Group-wise Operations
All techniques for this problem rely on the split-apply-combine strategy
First, take the data (or some object) and split it into smaller datasets on the basis of some variable
Dataset A
x
y
z
1
1
1
2
2
1
3
3
1
4
1
2
5
2
2
6
3
2
Datasets B and C (Dataset A split according to "z")
x
y
z
x
y
z
1
1
1
4
1
2
2
2
1
5
2
2
3
3
1
6
3
2
Second, apply some function to each one of the smaller datasets/objects
Example function: mean of variables "x" and "y"
Datasets B' and C'
mean(x)
mean(y)
z
mean(x)
mean(y)
z
2
2
1
5
2
2
Third, combine the results into a larger dataset/object
Datasets B' and C'
mean(x)
mean(y)
z
mean(x)
mean(y)
z
2
2
1
5
2
2
Dataset A'
mean(x)
mean(y)
z
2
2
1
5
2
2
Group-wise Operations/plyr
plyr is the go-to package for all your splitting-applying-combining needs
Among its many benefits (above base R capabilities):
Don't have to worry about different name, argument, or output consistencies
Easily parallelized
Input from, and output to, data frames, matricies, and lists
Progress bars for lengthy computation
Informative error messages
# Install the 'plyr' package (only necessary one time)# install.packages('plyr') # Not Run# Load the 'plyr' package (necessary every new R session)library(plyr)
Group-wise Operations/plyr/selecting functions
Two essential questions:
What is the class of your input object?
What is the class of your desired output object?
If you want to split a data frame, and return results as a data frame, you use ddply
If you want to split a data frame, and return results as a list, you use dlply
If you want to split a list, and return results as a data frame, you use ldply
Group-wise Operations/plyr/writing commands
All of the major plyr functions have the same basic syntax
xxply(.data = , .variables = , .fun = , ...)
Consider the case where we want to calculate vote choice statistics across race from a data.frame, and return them as a data.frame:
# Check for (and remove) factor levels in 'race' and 'pres04'is.factor(x = red.blue$race)
## $asian
## kerry bush nader other
## 1 61.67 38.13 1.459 0.9728
##
## $black
## kerry bush nader other
## 1 87.18 12.39 1.675 1.186
##
## $`hispanic/latino`
## kerry bush nader other
## 1 60.93 37.77 1.535 1.874
##
## $other
## kerry bush nader other
## 1 59.4 37.9 2.471 3.032
##
## $white
## kerry bush nader other
## 1 43.87 55.2 1.327 1.314
##
## $`NA`
## kerry bush nader other
## 1 55.07 45.09 2.92 3.001
##
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
## race
## 1 asian
## 2 black
## 3 hispanic/latino
## 4 other
## 5 white
## 6 <NA>
Consider the case where we want to calculate vote choice statistics across race from a list, and return them as a data.frame:
# Split the data.frame into a list on the basis of 'race'
red.blue.race.list <-split(x = red.blue, f = red.blue$race)
# Check the class of the objectclass(red.blue.race.list)
## [1] "list"
# Check for list element namesobjects(red.blue.race.list)
## $asian
## kerry bush nader other
## 1 61.67 38.13 1.459 0.9728
##
## $black
## kerry bush nader other
## 1 87.18 12.39 1.675 1.186
##
## $`hispanic/latino`
## kerry bush nader other
## 1 60.93 37.77 1.535 1.874
##
## $other
## kerry bush nader other
## 1 59.4 37.9 2.471 3.032
##
## $white
## kerry bush nader other
## 1 43.87 55.2 1.327 1.314
Group-wise Operations/plyr/functions
plyr can accomodate any user-defined function, but it also comes with some pre-defined functions that assist with the most common split-apply-combine tasks
We've already seen summarize, which creates user-specified vectors and combines them into a data.frame. Here are some other helpful functions:
transform: applies a function to a data.frame and adds new vectors (columns) to it
# Add a column containing the average age of the race of the individual
red.blue.transform <-ddply(.data = red.blue, .variables = .(race), .fun = transform,
race.avg.age =mean(x = age9, na.rm = T))
unique(red.blue.transform$race.avg.age)
## [1] 4.133 4.003 3.704 4.029 4.674 5.464
Note that transform can't do transformations that involve the results of other transformations from the same call
# Attempt to add new columns that draw on other (but still new) columns
red.blue.transform <-ddply(.data = red.blue, .variables = .(race), .fun = transform,
race.avg.age =mean(x = age9, na.rm = T), race.avg.age.plusone = race.avg.age +1)
## Error: object 'race.avg.age' not found
For this, we need mutate: just like transform, but it executes the commands iteratively so transformations can be carried out that rely on previous transformations from the same call
# Attempt to add new columns that draw on other (but still new) columns
red.blue.mutate <-ddply(.data = red.blue, .variables = .(race), .fun = mutate,
race.avg.age =mean(x = age9, na.rm = T), race.avg.age.plusone = race.avg.age +1)
unique(red.blue.mutate$race.avg.age)
## [1] 4.133 4.003 3.704 4.029 4.674 5.464
unique(red.blue.mutate$race.avg.age.plusone)
## [1] 5.133 5.003 4.704 5.029 5.674 6.464
Another very useful function is arrange, which orders a data frame on the basis of column contents
## state pres04 sex race age9
## 1 2 2 male white 1
## 2 2 2 female white 1
## 3 2 2 female white 1
## 4 2 2 male black 1
## 5 2 2 male white 1
## 6 2 1 male white 1
## 7 2 1 female hispanic/latino 1
## 8 2 1 male black 1
## 9 2 1 male white 1
## 10 2 1 male white 1
## 11 2 1 female white 1
## 12 2 2 male hispanic/latino 1
## 13 2 2 male white 1
## 14 2 1 male white 1
## 15 2 1 female white 1
## 16 2 1 female white 1
## 17 2 2 female white 1
## 18 2 2 male other 1
## 19 2 1 female hispanic/latino 1
## 20 2 NA female hispanic/latino 1
## 21 2 2 male white 1
## 22 2 1 female white 1
## 23 2 2 female white 1
## 24 2 1 female white 1
## 25 2 1 male white 1
## state pres04 sex race age9
## 1 2 9 male white 9
## 2 2 1 female white 9
## 3 2 2 male white 9
## 4 2 1 male white 9
## 5 2 2 female white 9
## 6 2 1 male white 9
## 7 2 2 female white 9
## 8 2 1 female white 9
## 9 2 2 female white 9
## 10 2 2 male white 9
## 11 2 2 female white 9
## 12 2 1 female white 9
## 13 2 2 female asian 9
## 14 2 2 male white 9
## 15 2 2 male white 9
## 16 2 1 female white 9
## 17 1 2 male white 9
## 18 1 1 male white 9
## 19 1 1 female <NA> 9
## 20 1 2 male white 9
## 21 1 2 female white 9
## 22 1 2 male white 9
## 23 1 1 male white 9
## 24 1 2 female white 9
## 25 1 NA female white 9
# Arrange by 'age9' then 'sex'
red.blue.age.sex <-plyr::arrange(df = red.blue, age9, sex)
red.blue.age.sex[1:25, 1:5]
## state pres04 sex race age9
## 1 2 2 male white 1
## 2 2 2 male black 1
## 3 2 2 male white 1
## 4 2 1 male white 1
## 5 2 1 male black 1
## 6 2 1 male white 1
## 7 2 1 male white 1
## 8 2 2 male hispanic/latino 1
## 9 2 2 male white 1
## 10 2 1 male white 1
## 11 2 2 male other 1
## 12 2 2 male white 1
## 13 2 1 male white 1
## 14 2 2 male white 1
## 15 2 2 male white 1
## 16 2 2 male white 1
## 17 2 2 male hispanic/latino 1
## 18 2 1 male white 1
## 19 2 2 male white 1
## 20 2 1 male white 1
## 21 2 2 male white 1
## 22 2 2 male white 1
## 23 2 2 male white 1
## 24 2 9 male white 1
## 25 2 2 male white 1
# Arrange by 'sex' (descending) then 'age9'
red.blue.sex.age <-plyr::arrange(df = red.blue, desc(sex), age9)
red.blue.sex.age[1:25, 1:5]
## state pres04 sex race age9
## 1 2 2 female white 1
## 2 2 2 female white 1
## 3 2 1 female hispanic/latino 1
## 4 2 1 female white 1
## 5 2 1 female white 1
## 6 2 1 female white 1
## 7 2 2 female white 1
## 8 2 1 female hispanic/latino 1
## 9 2 NA female hispanic/latino 1
## 10 2 1 female white 1
## 11 2 2 female white 1
## 12 2 1 female white 1
## 13 2 2 female white 1
## 14 2 1 female white 1
## 15 2 2 female white 1
## 16 2 2 female white 1
## 17 2 NA female white 1
## 18 2 1 female white 1
## 19 2 2 female other 1
## 20 2 1 female white 1
## 21 2 2 female white 1
## 22 2 2 female <NA> 1
## 23 2 2 female white 1
## 24 2 3 female white 1
## 25 2 9 female white 1
Group-wise Operations/dplyr
While plyr is a really great split-apply-combine tool, it can be a little slow when working over very large datasets. Consider this dataset containing all votes cast by all members of the U.S. House of Representatives since 1990 (roughly 7 million observations).
# CAUTION: this is a really big fileload(file ="datasets/ushouse_votes")
# View first six rows, check dimensionalityhead(ushouse.votes)
## name party state vote year vote.num yea.vote
## 1 Ackerman D NY Present 1990 1 0
## 2 Akaka D HI Present 1990 1 0
## 3 Alexander D AR Present 1990 1 0
## 4 Anderson D CA Present 1990 1 0
## 5 Andrews (TX) D TX Present 1990 1 0
## 6 Annunzio D IL Present 1990 1 0
dim(ushouse.votes)
## [1] 6754366 7
For extraordinarily large datasets, there is a new iteration of plyr being developed called dplyr. It is built for working with data.frames only and (currently) only accepts summarise(), mutate(), and arrange(), as well as select() for looking at certain variables and filter() for looking at certain rows.
# Install the 'dplyr' package (only necessary one time)# install.packages('dplyr') # Not Run# Load the 'dplyr' package (necessary every new R session)library(dplyr)
# Basic syntax for a dplyr call:# 1) Create a tbl (tabular data structure) that groups the data
data <-group_by(x = , ...)
# 2) Execute the pre-defined functionsummarise(.data = , ...)
A side-by-side comparison of ddply and dplyr (with benchmarks):
## state avg.yea
## 1 AK 0.5346
## 2 AL 0.5960
## 3 AR 0.6135
## 4 AZ 0.6048
## 5 CA 0.6000
## 6 CO 0.6084
## 7 CT 0.6322
## 8 DE 0.6375
## 9 FL 0.6040
## 10 GA 0.6121
## 11 HI 0.6091
## 12 IA 0.6213
## 13 ID 0.6057
## 14 IL 0.6009
## 15 IN 0.6117
## 16 KS 0.6254
## 17 KY 0.6160
## 18 LA 0.5890
## 19 MA 0.6127
## 20 MD 0.6274
## 21 ME 0.6342
## 22 MI 0.6111
## 23 MN 0.6135
## 24 MO 0.5985
## 25 MS 0.6116
## 26 MT 0.5923
## 27 NC 0.6222
## 28 ND 0.6412
## 29 NE 0.6233
## 30 NH 0.6544
## 31 NJ 0.6156
## 32 NM 0.6170
## 33 NV 0.6329
## 34 NY 0.6098
## 35 OH 0.6139
## 36 OK 0.6017
## 37 OR 0.6122
## 38 PA 0.6092
## 39 RI 0.6160
## 40 SC 0.6193
## 41 SD 0.6641
## 42 TN 0.6229
## 43 TX 0.5958
## 44 UT 0.6163
## 45 VA 0.6125
## 46 VT 0.6072
## 47 WA 0.6161
## 48 WI 0.6307
## 49 WV 0.6040
## 50 WY 0.5404
## 51 XX 0.3112
# dplyr
state.groups <-group_by(x = ushouse.votes, state)
time.dplyr <-system.time(avg.yea.states.dplyr <-summarise(.data = state.groups,
avg.yea =mean(x = yea.vote, na.rm = T)))
# Compare the execution times (ratio of ddply-to-plyr)
time.ddply/time.dplyr
## user system elapsed
## 16.50 18.01 18.31
Reshaping Data/reshape2
Often times, even before we're interested in doing all this group-wise stuff, we need to reshape our data. For instance, datasets often arrive at your desk in wide (long) form and you need to convert them to long (wide) form.
Though base R does have commands for reshaping data (including aggregate, by, tapply, etc.), each of their input commands are slightly different and are only suited for specific reshaping tasks.
The reshape2 package overcomes these argument and task inconsistencies to provide a simple, relatively fast way to alter the form of a data.frame. The package contains effectively two commands, and their functions are in their names: melt and cast
# Install the 'reshape2' package (only necessary one time)# install.packages('reshape2') # Not Run# Load the 'reshape2' package (necessary every new R session)library(reshape2)
Reshaping Data/reshape2/melt
melt() is used to convert wide-form data to long-form. The basic idea is to take your data.frame and melt it down to a minimal number of columns using two essential pieces of information:
Unit-of-Analysis identifiers, or columns you don't want to melt down
Characteristic variables, or columns you do want to melt down
## state pres04 sex race age9
## 1 2 1 female white 2
## 2 2 2 male white 1
## 3 2 1 female black 3
## 4 2 1 female black 3
## 5 2 1 female white 4
## 6 2 1 female white 3
## 7 2 1 male white 4
## 8 2 2 female white 1
## 9 2 2 female white 2
## 10 2 2 female white 1
Suppose we wanted to convert this data from its current wide format to an entirely long format. How to proceed?
First, select which columns you want to keep (i.e. not melt). In this case, I'm interested in having individual voters as my unit of analysis. Unfortunately, there is no column containing an individual identification number in this data, so I'll just add one as "individual":
## individual variable value
## 1 1 state 2
## 2 2 state 2
## 3 3 state 2
## 4 4 state 2
## 5 5 state 2
## 6 6 state 2
## 7 7 state 2
## 8 8 state 2
## 9 9 state 2
## 10 10 state 2
## 11 1 pres04 1
## 12 2 pres04 2
## 13 3 pres04 1
## 14 4 pres04 1
## 15 5 pres04 1
## 16 6 pres04 1
## 17 7 pres04 1
## 18 8 pres04 2
## 19 9 pres04 2
## 20 10 pres04 2
## 21 1 sex female
## 22 2 sex male
## 23 3 sex female
## 24 4 sex female
## 25 5 sex female
## 26 6 sex female
## 27 7 sex male
## 28 8 sex female
## 29 9 sex female
## 30 10 sex female
## 31 1 race white
## 32 2 race white
## 33 3 race black
## 34 4 race black
## 35 5 race white
## 36 6 race white
## 37 7 race white
## 38 8 race white
## 39 9 race white
## 40 10 race white
## 41 1 age9 2
## 42 2 age9 1
## 43 3 age9 3
## 44 4 age9 3
## 45 5 age9 4
## 46 6 age9 3
## 47 7 age9 4
## 48 8 age9 1
## 49 9 age9 2
## 50 10 age9 1
# If you want to melt ALL columns that aren't ID variables, you can also# omit the 'measure.vars' argumentmelt(data = red.blue.subset, id.vars ="individual")
## individual variable value
## 1 1 state 2
## 2 2 state 2
## 3 3 state 2
## 4 4 state 2
## 5 5 state 2
## 6 6 state 2
## 7 7 state 2
## 8 8 state 2
## 9 9 state 2
## 10 10 state 2
## 11 1 pres04 1
## 12 2 pres04 2
## 13 3 pres04 1
## 14 4 pres04 1
## 15 5 pres04 1
## 16 6 pres04 1
## 17 7 pres04 1
## 18 8 pres04 2
## 19 9 pres04 2
## 20 10 pres04 2
## 21 1 sex female
## 22 2 sex male
## 23 3 sex female
## 24 4 sex female
## 25 5 sex female
## 26 6 sex female
## 27 7 sex male
## 28 8 sex female
## 29 9 sex female
## 30 10 sex female
## 31 1 race white
## 32 2 race white
## 33 3 race black
## 34 4 race black
## 35 5 race white
## 36 6 race white
## 37 7 race white
## 38 8 race white
## 39 9 race white
## 40 10 race white
## 41 1 age9 2
## 42 2 age9 1
## 43 3 age9 3
## 44 4 age9 3
## 45 5 age9 4
## 46 6 age9 3
## 47 7 age9 4
## 48 8 age9 1
## 49 9 age9 2
## 50 10 age9 1
Note that melt collapses all of the measure variables into two columns: one containing the column/measurement name, the other containing the column/measurement value for that row. By default, these columns are named "variable" and "value", though they can be customized using the "variable.name" and "value.name" arguments. For example:
## individual characteristic response
## 1 1 state 2
## 2 2 state 2
## 3 3 state 2
## 4 4 state 2
## 5 5 state 2
## 6 6 state 2
## 7 7 state 2
## 8 8 state 2
## 9 9 state 2
## 10 10 state 2
## 11 1 pres04 1
## 12 2 pres04 2
## 13 3 pres04 1
## 14 4 pres04 1
## 15 5 pres04 1
## 16 6 pres04 1
## 17 7 pres04 1
## 18 8 pres04 2
## 19 9 pres04 2
## 20 10 pres04 2
## 21 1 sex female
## 22 2 sex male
## 23 3 sex female
## 24 4 sex female
## 25 5 sex female
## 26 6 sex female
## 27 7 sex male
## 28 8 sex female
## 29 9 sex female
## 30 10 sex female
## 31 1 race white
## 32 2 race white
## 33 3 race black
## 34 4 race black
## 35 5 race white
## 36 6 race white
## 37 7 race white
## 38 8 race white
## 39 9 race white
## 40 10 race white
## 41 1 age9 2
## 42 2 age9 1
## 43 3 age9 3
## 44 4 age9 3
## 45 5 age9 4
## 46 6 age9 3
## 47 7 age9 4
## 48 8 age9 1
## 49 9 age9 2
## 50 10 age9 1
Note also that one need not melt down all columns that aren't serving as ID columns. The melted data.frame will only contain the values of the measure variables you select. For instance:
## individual variable value
## 1 1 pres04 1
## 2 2 pres04 2
## 3 3 pres04 1
## 4 4 pres04 1
## 5 5 pres04 1
## 6 6 pres04 1
## 7 7 pres04 1
## 8 8 pres04 2
## 9 9 pres04 2
## 10 10 pres04 2
## 11 1 sex female
## 12 2 sex male
## 13 3 sex female
## 14 4 sex female
## 15 5 sex female
## 16 6 sex female
## 17 7 sex male
## 18 8 sex female
## 19 9 sex female
## 20 10 sex female
Reshaping Data/reshape2/cast
There are two main cast functions in the reshape2 package for converting data from a long format to a wide format: acast() (for producing arrays) and dcast() (for producing data frames)
## state pres04 sex race age9 individual
## 1 2 1 female white 2 1
## 2 2 2 male white 1 2
## 3 2 1 female black 3 3
## 4 2 1 female black 3 4
## 5 2 1 female white 4 5
## 6 2 1 female white 3 6
## 7 2 1 male white 4 7
## 8 2 2 female white 1 8
## 9 2 2 female white 2 9
## 10 2 2 female white 1 10
# Cast a data.frame containing the individual column and columns containing# the expansion of 'age9' on the basis of its unique valuesdcast(data = red.blue.subset, formula = individual ~age9, value.var ="age9")
## individual 1 2 3 4
## 1 1 NA 2 NA NA
## 2 2 1 NA NA NA
## 3 3 NA NA 3 NA
## 4 4 NA NA 3 NA
## 5 5 NA NA NA 4
## 6 6 NA NA 3 NA
## 7 7 NA NA NA 4
## 8 8 1 NA NA NA
## 9 9 NA 2 NA NA
## 10 10 1 NA NA NA
# Previously melted data
red.blue.melt
## individual variable value
## 1 1 pres04 1
## 2 2 pres04 2
## 3 3 pres04 1
## 4 4 pres04 1
## 5 5 pres04 1
## 6 6 pres04 1
## 7 7 pres04 1
## 8 8 pres04 2
## 9 9 pres04 2
## 10 10 pres04 2
## 11 1 sex female
## 12 2 sex male
## 13 3 sex female
## 14 4 sex female
## 15 5 sex female
## 16 6 sex female
## 17 7 sex male
## 18 8 sex female
## 19 9 sex female
## 20 10 sex female
# Cast a new data.frame from melted data.frame containing the individual# column and expanding the 'variable' columndcast(data = red.blue.melt, formula = individual ~variable, value.var ="value")
Once we've imported our data, summarized it, carried out group-wise operations, and perhaps reshaped it, we may also like to attempt causal inference.
This often requires doing the following:
Dealing with missing values
Carrying out Classical Hypothesis Tests
Estimating Regressions
Carryingout Regression Diagnostics
Inference/Missing Values
Having missing values can hinder the quality of your inferences for a variety of reasons. Luckily, statisticians have developed methods for dealing with missing values by imputing missing values. One of the best packages for this is Amelia.
# Install the 'Amelia' package (only necessary one time)# install.packages('Amelia') # Not Run# Load the 'Amelia' package (necessary every new R session)library(Amelia)
Consider the following panel dataset containing economic and demographic data for several African countries over roughly two decades starting in the 1970s.
data(africa) # Load the 'africa' data.set (included in Amelia)# Check dimensionality, colnames, head, tail, and structuredim(africa)
# Most important object is likely 'imputations' (this contains a list of# data.frames with imputed values included). Must be sure to apply your# analyses (e.g. regressions) over each data.frame in the list, and then# aggregate the resultsstr(africa.imputed$imputations)
# Can also apply summary methods to objects of class 'amelia'summary(africa.imputed)
##
## Amelia output with 5 imputed datasets.
## Return code: 1
## Message: Normal EM convergence.
##
## Chain Lengths:
## --------------
## Imputation 1: 3
## Imputation 2: 2
## Imputation 3: 3
## Imputation 4: 3
## Imputation 5: 2
##
## Rows after Listwise Deletion: 115
## Rows after Imputation: 120
## Patterns of missingness in the data: 3
##
## Fraction Missing for original variables:
## -----------------------------------------
##
## Fraction Missing
## year 0.00000
## country 0.00000
## gdp_pc 0.01667
## infl 0.00000
## trade 0.04167
## civlib 0.00000
## population 0.00000
An example of carring out regressions over each of the imputed datasets:
# Extract the imputed datasets from the amelia object
datasets <-africa.imputed$imputations
# Apply the regression over each element of the listlapply(X = datasets, FUN = function(x) {
summary(lm(formula = gdp_pc ~infl +trade, data = x))$coefficients
})
Suppose we have two different samples, A and B, both drawn from the standard normal distribution:
a <-rnorm(n =5000, mean =0, sd =1)
b <-rnorm(n =5000, mean =0, sd =1)
Suppose we also have a third sample, C, drawn from the normal distribution with mean=1 and sd=0:
c <-rnorm(n =5000, mean =1, sd =1)
One can test for differences in these distributions in either a) their means using t-tests, or b) their entire distributions using ks-tests
# Basic Callt.test(x = , y = , var.equal = , conf.level = , formula = )
ks.test(x = , y = )
t.test(x = a, y = b) # No difference in means
##
## Welch Two Sample t-test
##
## data: a and b
## t = 1.295, df = 9997, p-value = 0.1953
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.01315 0.06434
## sample estimates:
## mean of x mean of y
## 0.019855 -0.005742
t.test(x = a, y = c) # Difference in means
##
## Welch Two Sample t-test
##
## data: a and c
## t = -49.69, df = 9997, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.0299 -0.9517
## sample estimates:
## mean of x mean of y
## 0.01985 1.01067
ks.test(x = a, y = b) # No difference in distributions
##
## Two-sample Kolmogorov-Smirnov test
##
## data: a and b
## D = 0.0186, p-value = 0.3527
## alternative hypothesis: two-sided
ks.test(x = a, y = c) # Difference in distributions
##
## Two-sample Kolmogorov-Smirnov test
##
## data: a and c
## D = 0.3852, p-value < 2.2e-16
## alternative hypothesis: two-sided
Inference/Regression
Running regressions in R is extremely simple, very straightforwd (though doing things with standard errors requires a little extra work)
Most basic, catch-all regression function in R is glm
glm fits a generalized linear model with your choice of family/link function (gaussian, logit, poisson, etc.)
lm is just a standard linear regression (equivalent to glm with family=gaussian(link="identity"))
The basic glm call looks something like this:
glm(formula = y ~x1 +x2 +x3 +..., family =familyname(link ="linkname"),
data = )
There are a bunch of families and links to use (help(family) for a full list), but some essentials are binomial(link = "logit"), gaussian(link = "identity"), and poisson(link = "log")
Example: suppose we want to regress being an old man on political party identification, income, and religion using a logit model. The glm call would be something like this:
# Create an indicator for men aged 65 or over
red.blue$oldman <-ifelse(test = (red.blue$age65 == "65 or over" &red.blue$sex == "male"), yes =1, no =0)
# Regress being a 65+ year old male on partyid, income, and religion (logit# model)
oldman.reg <-glm(formula = oldman ~partyid +income +relign8, family =binomial(link ="logit"),
data = red.blue)
When we store this regression in an object, we get access to several items of interest
# View objects contained in the regression outputobjects(oldman.reg)
Note that, in our results, R has broken up our variables into their different factor levels (as it will do whenever your regressors have factor levels)
If your data aren't factorized, you can tell glm to factorize a variable (i.e. create dummy variables on the fly) by writing
glm(formula = y ~x1 +x2 +factor(x3), family =family(link ="link"), data = )
There are also some useful shortcuts for regressing on interaction terms:
x1:x2 interacts all terms in x1 with all terms in x2
summary(glm(formula = oldman ~partyid:income, family =binomial(link ="logit"),
data = red.blue))
Sometimes, there may be considerable uncertainty as to the proper form of the model to be estimated. For this, one can source the glm.sim function from Github using the package devtools.
# Install the 'devtools' package (only necessary one time)# install.packages('devtools') # Not Run# Load the 'devtools' package (necessary every new R session)library(devtools)
# Source the 'glm.sim' function from gist 9359460source_gist("https://gist.github.com/ckrogs/9359460")
Consider an example with the "africa" dataset. Suppose we wanted to regress "gdp_pc" on some combination of "infl", "trade", "civlib", and "population", but we were unsure which model best explained the outcome. glm.sim can runs all the possible models and aggregates the results.
# Remove extraneous columns from the data
africa.subset <-africa[, c(3:7)]
# Run glm.sim with graphs
output <-glm.sim(data = africa.subset, dv ="gdp_pc", graphs = T)
# Run glm.sim with graphs, best model, and robust standard errors
output <-glm.sim(data = africa.subset, dv ="gdp_pc", graphs = T, graph.best =1,
robust = T)
# Run glm.sim with graphs, best model, and first-order interactions
output <-glm.sim(data = africa.subset, dv ="gdp_pc", graphs = T, graph.best =1,
interactions = T)
# Run glm.sim with graphs, best model, and quadratic terms
output <-glm.sim(data = africa.subset, dv ="gdp_pc", graphs = T, graph.best =1,
quadratics = T)
3 main facilities for producing graphics in R: base, lattice, and ggplot2
In practice, these facilities are grouped into two camps: "basic" and "advanced"
A better formulation: quick/dirty v. involved/fancy
quick/dirty v. involved/fancy
Recall that R is an object-oriented programming language
tips <-reshape2::tips # Load dataset on tipping behavior included with reshape2 packageattributes(tips) # Check attributes of the tips dataset (names, row.names, class)
# Create an object of class 'lm' (linear model), regressing tip on some# covariates
tips.reg <-lm(formula = tip ~total_bill +sex +smoker +day +time +size,
data = tips)
attributes(tips.reg) # Check attributes of the tips.reg object (names, class)
# There are many variables, so it can be helpful to extract their classes# (mostly to check for factors) via a quick for-loop
classes <-NULL# Create a placeholder for the class output# for-loop
for (i in 1:ncol(data)) {
classes <-c(classes, class(data[, i]))
}
sort(classes) # Sort the results alphabetically (appears to be some factors, a few integers, and many numerics)
Note that data-bearing arguments of length creater than 1 ("x" in this case) must be vectorized in some manner (usually using the '$' column vector method)
plot(x = vturn) # Attempt to plot variable 'vturn' without vectorization
## Error: object 'vturn' not found
plot(x = data$vturn) # Attempt to plot variable 'vturn' with '$' vectorization
Basic call takes the following form
plot(x = , y = , type ="")
base/type (scatter, line, both)
The "type" argument accepts the following character indicators
plot(x = data$year, y = data$realgdpgr, type ="p") # Plot Year on x-axis and Real GDP Growth on y-axis
"l" -- line graphs
plot(x = data$year, y = data$realgdpgr, type ="l") # Plot Year on x-axis and Real GDP Growth on y-axis
"b" -- both line and point plots
plot(x = data$year, y = data$realgdpgr, type ="b") # Plot Year on x-axis and Real GDP Growth on y-axis
base/type (histograms, density plots)
Certain plot types require different calls outside of the "type" argument
Ex) Histograms
hist(x = data$vturn) # Plot histogram of voter turnout
hist(x = data$vturn, breaks =2) # Plot histogram of voter turnout, with 2 breaks
hist(x = data$vturn, breaks =50) # Plot histogram of voter turnout, with 50 breaks
Ex) Density plots
vturn.density <-density(x = data$vturn, na.rm = T) # Create a density object (NOTE: be sure to remove missing values)class(vturn.density) # Check the class of the object
## [1] "density"
vturn.density # View the contents of the object
##
## Call:
## density.default(x = data$vturn, na.rm = T)
##
## Data: data$vturn (1134 obs.); Bandwidth 'bw' = 2.826
##
## x y
## Min. : 26.5 Min. :0.00001
## 1st Qu.: 46.0 1st Qu.:0.00248
## Median : 65.4 Median :0.00626
## Mean : 65.4 Mean :0.01285
## 3rd Qu.: 84.8 3rd Qu.:0.02472
## Max. :104.3 Max. :0.03609
plot(x = vturn.density) # Plot the density object
plot(x =density(x = data$vturn, bw =2, na.rm = T)) # Plot the density object, bandwidth of 2
plot(x =density(x = data$vturn, bw =0.5, na.rm = T)) # Plot the density object, bandwidth of 0.5
plot(x =density(x = data$vturn, bw =6, na.rm = T)) # Plot the density object, bandwidth of 6
base/options [labeling]
Basic call with popular labeling arguments
plot(x = , y = , type ="", xlab ="", ylab ="", main ="")
From the previous example
plot(x = data$year, y = data$realgdpgr, type ="p", xlab ="Year", ylab ="Real GDP Growth",
main ="This Graph is Great") # Add labels for axes and overall plot
base/options [axis + size scaling]
Basic call with popular scaling arguments
plot(x = , y = , type ="", xlim = , ylim = , cex = )
From the previous example
plot(x = data$year, y = data$realgdpgr, type ="p") # Create a basic plot
# Limit the years (x-axis) to between 1976 and 1991plot(x = data$year, y = data$realgdpgr, type ="p", xlim =c(1976, 1991))
# Limit the years (x-axis) to between 1976 and 1991, increase point size to# 2plot(x = data$year, y = data$realgdpgr, type ="p", xlim =c(1976, 1991), cex =2)
# Limit the years (x-axis) to between 1976 and 1991, decrease point size to# 0.5plot(x = data$year, y = data$realgdpgr, type ="p", xlim =c(1976, 1991), cex =0.5)
base/options [graphical parameters]
Basic call with popular scaling arguments
plot(x = , y = , type ="", col ="", pch = , lty = , lwd = )
plot(x = data$year, y = data$realgdpgr, type ="p", pch =3) # Change point style to crosses
plot(x = data$year, y = data$realgdpgr, type ="p", pch =15) # Change point style to filled squares
# Change point style to filled squares and increase point size to 3plot(x = data$year, y = data$realgdpgr, type ="p", pch =15, cex =3)
plot(x = data$year, y = data$realgdpgr, type ="p", pch ="w") # Change point style to 'w'
# Change point style to '$' and increase point size to 2plot(x = data$year, y = data$realgdpgr, type ="p", pch ="$", cex =2)
Line Styles and Widths
# Create a data.frame containing yearly average Real GDP Growth over all# countrieslibrary(plyr)
# Split-apply-combine via plyr
mean.growth <-ddply(.data = data, .variables = .(year), summarize, mean =mean(realgdpgr))
head(mean.growth)
## year mean
## 1 1960 NA
## 2 1961 5.612
## 3 1962 4.830
## 4 1963 5.310
## 5 1964 6.618
## 6 1965 5.010
# Line plot with solid lineplot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =1)
# Line plot with medium dashed lineplot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =2)
# Line plot with short dashed lineplot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =3)
# Change line width to 2plot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =3, lwd =2)
# Change line width to 3plot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =3, lwd =3)
# Change line width to 4plot(x = mean.growth$year, y = mean.growth$mean, type ="l", lty =3, lwd =4)
base/layering
Layering is accomplished by plotting succesive commands of "lines()", "points()", etc. after "plot()"
# Subset data to create a few data.frames containing data for a single# country
france.growth <-data[data$country == "France", ]
italy.growth <-data[data$country == "Italy", ]
spain.growth <-data[data$country == "Spain", ]
Successive calls to "plot()" returns two different plots
# First call to plotplot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
# Second call to plotplot(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
But using "lines()" for the second and subsequent calls layers the plots
# First call to plotplot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
# First call to lineslines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
# Second call to lineslines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
The same is true for "points()"
# First call to plotplot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
# First call to lineslines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
# Second call to lineslines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
# First call to pointspoints(x =1986, y =6, pch =13, col =colors()[116])
# Second call to pointspoints(x =1986, y = -1, pch =13, col =colors()[24])
# Third call to pointspoints(x =1966, y =2, pch =13, col =colors()[391])
# Fourth call to pointspoints(x =2008, y =4, pch =13, col =colors()[8])
plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
# First call to texttext(x =1967, y = -1, labels ="No dot here")
Reference Lines
plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
text(x =1967, y = -1, labels ="No dot here")
# First call to ablineabline(v =1980, h =0)
Legends
plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
text(x =1967, y = -1, labels ="No dot here")
abline(v =1980, h =0)
# First call to legend (note the vector position correspondence between each# of the argument values)legend("topright", legend =c("France", "Italy", "Spain"), col =c("red", "blue",
"darkgreen"), lty =c(1, 1, 1), lwd =c(2, 2, 2), cex =0.8)
base/tables
Can form tables of graphs using the "par()" call like so:
par(mrow =c(ncols, nrows))
# STEP 1: Call 'par() for a 2x2 table'par(mfrow =c(2, 2))
# STEP 2: Plot #1plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
# STEP 3: Plot #2plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
text(x =1967, y = -1, labels ="No dot here")
# STEP 4: Plot #3plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
text(x =1967, y = -1, labels ="No dot here")
abline(v =1980, h =0)
# STEP 5: Plot #4plot(x = france.growth$year, y = france.growth$realgdpgr, type ="l", col ="red",
lwd =2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type ="l", col ="blue",
lwd =2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type ="l", col ="darkgreen",
lwd =2)
points(x =1986, y =6, pch =13, col =colors()[116])
points(x =1986, y = -1, pch =13, col =colors()[24])
points(x =1966, y =2, pch =13, col =colors()[391])
points(x =2008, y =4, pch =13, col =colors()[8])
text(x =1967, y = -1, labels ="No dot here")
abline(v =1980, h =0)
legend("topright", legend =c("France", "Italy", "Spain"), col =c("red", "blue",
"darkgreen"), lty =c(1, 1, 1), lwd =c(2, 2, 2), cex =0.8)
Out with the old...
Base graphics are really great, but they're not like this (created with a single line of code)
lattice (Deepayan Sarkar, ISI, Delhi)
ggplot2 (Hadley Wickham, again)
Both are built on "grid", both are really huge improvements over base R graphics
Both also have entire books written about them (~200-300 pp. each)
lattice v. ggplot2
lattice is
faster (though only noticeable over many and large plots)
simpler (at first)
better at trellis graphs
able to do 3d graphs
ggplot2 is
generally more elegant
more syntactically logical (and therefore simpler, once you learn it)
better at grouping
able to interface with maps
Basic usage: lattice
The general call for lattice graphics looks something like this:
graph_type(formula, data=, [options])
The specifics of the formula differ for each graph type, but the general format is straightforward
y # Show the distribution of y
y ~x # Show the relationship between x and y
y ~x |A # Show the relationship between x and y conditional on the values of A
y ~x |A *B # Show the relationship between x and y conditional on the combinations of A and B
z ~y *x # Show the 3D relationship between x, y, and z
Basic usage: ggplot2
The general call for ggplot2 graphics looks something like this:
# Create data.frame of average growth rates by country over time
growth <-ddply(.data = data, .variables = .(country), summarize, mean =mean(realgdpgr,
na.rm = T))
barchart(mean ~country, data = growth) # lattice
ggplot(data = volcano3d, aes(x = xvar, y = yvar, z = zvar)) +geom_contour() # ggplot2
lattice v. ggplot2: tile/image/level plots
levelplot(zvar ~xvar +yvar, data = volcano3d) # lattice
ggplot(data = volcano3d, aes(x = xvar, y = yvar, z = zvar)) +geom_tile(aes(fill = zvar)) # ggplot2
lattice: 3D plots
# Create a subset of the dataset containing only data for France
france.data <-data[data$country == "France", ]
cloud(outlays ~year *vturn, data = france.data)
# Create a subset of the dataset containing only data for Greece, Portugal,# Ireland, and Spain
pigs.data <-data[data$country %in%c("Greece", "Portugal", "Ireland", "Spain"),
]
cloud(outlays ~year *vturn |country, data = pigs.data)
ggplot2: Panel plots
ggplot(data = pigs.data, aes(x = year, y = vturn, color = country)) +geom_line()
lattice v. ggplot2: options [labeling]
xyplot(outlays ~vturn, data = data, xlab ="Voter Turnout (%)", ylab ="Government Outlays",
main ="This Graph is Also Great") # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) +geom_point() +xlab(label ="Voter Turnout (%)") +ylab(label ="Government Outlays") +ggtitle(label ="This Graph is Also Also Great") # ggplot2
lattice v. ggplot2: options [axis + size scaling]
xyplot(outlays ~vturn, data = data, xlim =c(80, 100)) # lattice
By now, you might be noticing some trends in how these two packages approach graphics
lattice tends to focus on a particular type of graph and how to represent cross-sectional variation by splitting it up into smaller chunks
Becoming a proficient user of lattice requires learning a huge array of graph-specific formulas and options
ggplot2 tries to represent much more of the cross-sectional variation by making use of various "aesthetics"; general approach is based on The Grammar of Graphics
ggplot2 and the Grammar of Graphics
Basic idea is that the visualization of all data requires four items
One or more statistics conveying information about the data (identities, means, medians, etc.)
A coordinate system that differentiates between the intersections of statistics (at most two for ggplot, three for lattice)
Geometries that differentiate between off-coordinate variation in kind
Scales that differentiate between off-coordinate variation in degree
ggplot2 allows the user to manipulate all four of these items
ggplot2 is optimized for showing variation on all four aesthetic types
# Differences in kind using colorggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_line(aes(color = country))
Note what happens when we specify the color parameter outside of the aesthetic operator. ggplot2 views these specifications as invalid graphical parameters.
# Differences in kind using line typesggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_line(aes(linetype = country))
# Differences in kind using point shapesggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point(aes(shape = country))
# Differences in degree using colorggplot(data = pigs.data, aes(x = year, y = realgdpgr)) +geom_point(aes(color = vturn))
# Differences in degree using point sizeggplot(data = pigs.data, aes(x = year, y = realgdpgr)) +geom_point(aes(size = vturn))
# Multiple non-cartesian aesthetics (differences in kind using color, degree# using point size)ggplot(data = pigs.data, aes(x = year, y = realgdpgr)) +geom_point(aes(color = country,
size = vturn))
Fitted lines and curves with ggplot2
ggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point()
# Add linear model (lm) smootherggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point() +geom_smooth(method ="lm")
# Add local linear model (loess) smoother, span of 0.75ggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point() +geom_smooth(method ="loess",
span =0.75)
# Add local linear model (loess) smoother, span of 0.25ggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point() +geom_smooth(method ="loess",
span =0.25)
# Add linear model (lm) smoother, no standard error shadingggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point() +geom_smooth(method ="lm",
se = F)
# Add local linear model (loess) smoother, no standard error shadingggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point() +geom_smooth(method ="loess",
se = F)
# Add a local linear (loess) smoother for each countryggplot(data = pigs.data, aes(x = year, y = vturn)) +geom_point(aes(color = country)) +geom_smooth(aes(color = country))
# Add a local linear (loess) smoother for each country, no standard error# shadingggplot(data = pigs.data, aes(x = year, y = realgdpgr)) +geom_point(aes(color = country,
size = vturn)) +geom_smooth(aes(color = country), se = F)
lattice v. ggplot2: tables
Both lattice and ggplot2 graphs can be combined using the grid.arrange() function in the gridExtra package
# Initialize gridExtra librarylibrary(gridExtra)
# Create 3 plots to combine in a table
plot1 <-ggplot(data = pigs.data, aes(x = year, y = vturn, color = )) +geom_line(aes(color = country))
plot2 <-ggplot(data = pigs.data, aes(x = year, y = vturn, linetype = )) +geom_line(aes(linetype = country))
plot3 <-ggplot(data = pigs.data, aes(x = year, y = vturn, shape = )) +geom_point(aes(shape = country))
# Call grid.arrangegrid.arrange(plot1, plot2, plot3, nrow =3, ncol =1)
Exporting
Two basic image types
Raster/Bitmap (.png, .jpeg)
Every pixel of a plot contains its own separate coding; not so great if you want to resize the image
